


### Figure 4. Univariate LM test statistics across varying sample sizes


## Package required 

install.packages("lavaan")
install.packages("tidyverse")

library(lavaan)
library(tidyverse)


###### Population model setup 



popmodel<-'
f1 =~ 0.65*x1 + 0.65*x2 + 0.7*x3 + 0.7*x4 + 0.7*x5 + 0.7*x6 + 0.6*x7 + 0.5*x8 + 0.5*x9
f2 =~ 0.6*x9 + 0.6*x10 + 0.6*x11 + 0.7*x12 + 0.7*x13 + 0.5*x14 + 0.5*x15 + 0.65*x16 + 0.5*x6 + 0.5*x20
f3 =~ 0.5*x17 + 0.5*x18 + 0.5*x19 + 0.6*x20 + 0.6*x21 + 0.6*x22 + 0.7*x23 + 0.7*x24 + 0.45*x16

f3 ~~ 0.4*f1 + 0.5*f2
f1 ~~ 0.3*f2

f1 ~ 0*1
f2 ~ 0*1
f3 ~ 0*1

x1 ~~ 0.5775*x1
x2 ~~ 0.5775*x2

x3 ~~ 0.51*x3
x4 ~~ 0.51*x4
x5 ~~ 0.51*x5
x6 ~~ 0.05*x6

x7 ~~ 0.64*x7
x8 ~~ 0.21*x8
x9 ~~ 0.21*x9

x10 ~~  0.64*x10
x11 ~~  0.51*x11

x12 ~~ 0.51*x12
x13 ~~ 0.75*x13
x14 ~~ 0.75*x14
x15 ~~ 0.5775*x15

x16 ~~ 0.7975*x16
x17 ~~ 0.75*x17

x18 ~~ 0.75*x18
x19 ~~ 0.25*x19
x20 ~~ 0.64*x20
x21 ~~ 0.64*x21
x22 ~~ 0.64*x22

x23 ~~ 0.51*x23
x24 ~~ 0.51*x24


'


#### Analysis model setup

mis_model<-'  

f1 =~ x1 + x2 + x3 + x4 + x5 + x6  + x7 + x8
f2 =~ x9 + x10 + x11 + x12 + x13 + x14 + x15 + x16
f3 =~ x17 + x18 + x19 + x20 + x21 + x22 + x23 + x24
f3 ~~ f1
f1 ~~ f2


'

###  Test the model in different sample sizes

##################################################

                 N=100

##################################################

bin_1000 <- data.frame(id = character(), mi = numeric(), trial = integer())

#seed number is the same as N
set.seed(103)

for (i in 1:500) {
  
  mydata <- simulateData(popmodel, model.type = "sem", sample.nobs = 100L) # simulate data for each trial
  
  fit1 <- sem(model = mis_model, data = mydata, meanstructure = FALSE, likelihood = "wishart", estimator = "ML")
  
  mi_df <- modindices(fit1, minimum.value = 0, sort = TRUE) %>%
    
    rownames_to_column(var = "id") # add an ID column to mi_df
  
  top20 <- mi_df[1:20, c("id", "mi")]
  
  top20$mi <- round(top20$mi, 3) # round mi to 3 decimal places
  
  top20$trial <- i
  
  bin_1000 <- rbind(bin_1000, top20)
  
}



bin_1000



par(mar = c(5, 4, 4, 4) + 0.1)
highlighted_bars <- c(76, 55, 82, 102)
colors <- ifelse(as.numeric(bin_1000$id) %in% highlighted_bars, "blue", "gray")
plot(as.numeric(bin_1000$id), as.numeric(bin_1000$mi), xaxt = 'n', xlab = "Parameters", ylab = "Univariate LM Test", col = colors, main = "N = 100")



##################################################

                N=500

##################################################


bin_1000 <- data.frame(id = character(), mi = numeric(), trial = integer())

set.seed(500)

for (i in 1:500) {
  
  mydata <- simulateData(popmodel, model.type = "sem", sample.nobs = 500L) # simulate data for each trial
  
  fit1 <- sem(model = mis_model, data = mydata, meanstructure = FALSE, likelihood = "wishart", estimator = "ML")
  
  mi_df <- modindices(fit1, minimum.value = 0, sort = TRUE) %>%
    
    rownames_to_column(var = "id") # add an ID column to mi_df
  
  top20 <- mi_df[1:20, c("id", "mi")]
  
  top20$mi <- round(top20$mi, 3) # round mi to 3 decimal places
  
  top20$trial <- i
  
  bin_1000 <- rbind(bin_1000, top20)
  
}



bin_1000



par(mar = c(5, 4, 4, 4) + 0.1)
highlighted_bars <- c(76, 55, 82, 102)
colors <- ifelse(as.numeric(bin_1000$id) %in% highlighted_bars, "blue", "gray")
plot(as.numeric(bin_1000$id), as.numeric(bin_1000$mi), xaxt = 'n', xlab = "Parameters", ylab = "Univariate LM Test", col = colors, main = "N = 500")






##################################################

N=1000

##################################################

library(tidyverse)


bin_1000 <- data.frame(id = character(), mi = numeric(), trial = integer())

#seed number is the same as N
set.seed(1000)

for (i in 1:500) {
  
  mydata <- simulateData(popmodel, model.type = "sem", sample.nobs = 1000L) # simulate data for each trial
  
  fit1 <- sem(model = mis_model, data = mydata, meanstructure = FALSE, likelihood = "wishart", estimator = "ML")
  
  mi_df <- modindices(fit1, minimum.value = 0, sort = TRUE) %>%
    
    rownames_to_column(var = "id") # add an ID column to mi_df
  
  top20 <- mi_df[1:20, c("id", "mi")]
  
  top20$mi <- round(top20$mi, 3) # round mi to 3 decimal places
  
  top20$trial <- i
  
  bin_1000 <- rbind(bin_1000, top20)
  
}



bin_1000

par(mar = c(5, 4, 4, 4) + 0.1)
highlighted_bars <- c(76, 55, 82, 102)
colors <- ifelse(as.numeric(bin_1000$id) %in% highlighted_bars, "blue", "gray")
plot(as.numeric(bin_1000$id), as.numeric(bin_1000$mi), xaxt = 'n', xlab = "Parameters", ylab = "Univariate LM Test", col = colors, main = "N = 1000")



##################################################

              N=5000

##################################################

bin_1000 <- data.frame(id = character(), mi = numeric(), trial = integer())

set.seed(5000)

for (i in 1:500) {
  
  mydata <- simulateData(popmodel, model.type = "sem", sample.nobs = 5000L) # simulate data for each trial
  
  fit1 <- sem(model = mis_model, data = mydata, meanstructure = FALSE, likelihood = "wishart", estimator = "ML")
  
  mi_df <- modindices(fit1, minimum.value = 0, sort = TRUE) %>%
    
    rownames_to_column(var = "id") # add an ID column to mi_df
  
  top20 <- mi_df[1:20, c("id", "mi")]
  
  top20$mi <- round(top20$mi, 3) # round mi to 3 decimal places
  
  top20$trial <- i
  
  bin_1000 <- rbind(bin_1000, top20)
  
}


bin_1000

par(mar = c(5, 4, 4, 4) + 0.1)
highlighted_bars <- c(76, 55, 82, 102)
colors <- ifelse(as.numeric(bin_1000$id) %in% highlighted_bars, "blue", "gray")
plot(as.numeric(bin_1000$id), as.numeric(bin_1000$mi), xaxt = 'n', xlab = "Parameters", ylab = "Univariate LM Test", col = colors, main = "N = 5000")



### End of code
